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Abstract 

R.  M.  Lewitt  has  introduced  a family  of  compactly  supported  radial  basis  functions 
which  are  particularly  useful  in  discretising  for  inversion  ill-posed  problems  involving  line 
integrals.  We  consider  some  practical  considerations  in  their  use  and  implementation, 
compare  square  and  triangular  grids  of  the  functions  in  two  dimensions,  and  describe 
some  particularly  favourable  choices  of  the  defining  parameters. 

1 Introduction 

In  the  article  [5],  R.  M.  Lewitt  introduced  a family  of  window  functions 

Mr)  = / (1  - - (r/a)2)^)/Im(a),  0 < r < a, 

• w \ 0,  r>  a,  K ' ’ 

where  Im  is  the  modified  Bessel  function  of  order  m (see  Ch.  Ill,  3.7  [13]).  The  implicit 
dependence  of  ip  on  the  parameters  a > 0,  a > 0 and  m € N is  discussed  below.  Lewitt ’s 
motivation  for  studying  these  functions  is  the  use  of  translates  of  the  radially  symmetric 
function 

W(x)  = tHIMI)  (x  G Rrf) 

(see  Figure  1)  as  a basis  for  the  discretisation  of  tomographic  problems  [8,  9].  Such  a basis 
overcomes  a number  of  difficulties  associated  with  the  usual,  pixel-based,  representation 
in  problems  involving  the  recovery  of  function  from  a set  of  line,  curve  or  strip  integrals 
across  its  domain,  while  retaining  the  advantage  of  a sparse  discretisation.  The  author’s 
interest  in  these  functions  arises  in  their  application  to  a Radon-like  problem  in  the 
remote  sensing  of  ocean  waves  [15],  a detailed  exposition  of  which  may  be  found  in  [3]. 

2 Discretising  x-ray  problems 

The  discretisation  of  an  x-ray  transform  inversion  problem  with  Lewitt ’s  basis  is  straight- 
forward. Given  a set  of  centres  x*  € Rd,  one  represents  the  (unknown)  function  / as  a 
linear  combination  of  the  translates  of  T, 

f(x)  = '^/^(x-Xi)  (x  G Rd).  (2.1) 
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Fig.  1.  Lewitt’s  radial  basis  function  in  dimension  2 with  m — 2,  a = 3. 

The  given  data  in  such  problems  are  the  values  Ij  of  integrals  of  / over  lines  (or  more 
generally,  submanifolds)  Lj 

Ij  - [ f(x)  = 

JLj  i 

The  latter  integral  in  (2.2)  is  the  projection  or  Abel  transform  of  /,  which  can  be  calcu- 
lated explicitly  in  the  linear  case.  For  a line  Lj  whose  closest  point  to  xt  is  at  a distance 
s from  it,  and  with  the  dependence  of  if  on  m here  made  explicit, 

2/  i’mW s2  + t2)dt  ==aIm+1/2^. 

JO  \ & ) 

(see  A7,  [5]).  Thus  (2.2)  reduces  to  a linear  system  which  may  be  solved  for  the  coefficients 
&.  If  the  support  of  the  basis  functions  is  small  (i.  e.,  if  a is  small)  then  this  linear  system 
has  an  unstructured  sparsity  which  can  be  exploited  by,  for  example,  an  iterative  row- 
action  solution  method  [2]. 

The  computational  cost  of  such  a discretisation  lies  mainly  in  the  evaluation  of  the 
Abel  transform  which  requires  the  calculation  of  a Bessel  function.  Fortunately,  Bessel 
functions  of  half-integer  order  can  be  calculated  efficiently  from  their  recurrence  relations 
(see  the  Atlas,  [12],  for  details). 

The  discretisation  techniques  describe  here  can  also  be  applied  to  problems  in  which 
the  integrals  are  over  curves  of  sufficient  smoothness  to  allow  a local  linear  approxima- 
tion. 
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3 Fourier  transform  and  invertibility 

The  Fourier  transform  of  the  (/-dimensional  basis  function  T„,  is  radially  symmetric  and 
given  in  (A3)  of  [5]  as 


fmW  = 


adam( 2ir)d/2  Jm+d/2{z) 
Im(a)  zm+d/2  ’ 


2 = y/  (27t«||.t||)2  — a2. 


(3.1) 


The  presence  of  the  Bessel  function  Jm+d/ 2(2)  in  this  expression  clearly  implies  that  it  is 
not  non-negative,  and  so  by  Bochner’s  characterisation  of  positive  definite  translation- 
invariant  functions,  $ is  not  positive  definite  for  any  choices  of  the  parameters. 

This  fact  denies  us  the  attractive  approximation  theory  of  the  compactly  supported 
radial  functions  of  Wu,  Wendland  and  Buhmann  (Section  3,  [1]).  In  particular,  there  is 
no  guarantee,  per  se,  on  the  invertibility  of  the  interpolation  matrix  [T(.t,'  — Xj)\,  needed 
to  ensure  that  (2.1)  can  represent  an  arbitrary  function  at  its  centres.  However,  this 
interpolation  matrix  is  invertible  if  it  is  strictly  diagonally  dominant  (Corollary  5.6.17, 
[4])  which,  for  a set  of  centres  on  a uniform  grid  T,  holds  if 


*(0)>  E *(*)•  (3-2) 

*er\{o} 

Values  of  the  parameters  for  which  (3.2)  is  satisfied  for  the  square  planar  grid  AZ2  are 
shown  in  Figure  2. 

As  is  noted  in  [5] , there  are  several  reasons  why  a rapid  decay  of  the  Fourier  transform 
of  the  basis  function  is  advantageous  in  functional  representation  for  the  inversion  of  x- 
ray  and  related  transforms. 

• Such  inversions  may  be  complicated  by  functions  in  the  nullspace  of  the  transform, 
so-called  ghosts.  For  some  transforms  [7]  it  can  be  shown  that  such  functions  have 
a Fourier  transform  which  is  small  close  to  zero  in  the  frequency  domain,  and  so 
representation  by  a basis  with  Fourier  transform  localised  around  zero  will  suppress 
these  ghosts. 

• These  inversions  are  often  ill-posed  and  the  given  data  noisy.  Representation  of  the 
sought  function  by  a basis  with  localised  Fourier  transform  imposes  smoothness, 
and  so  acts  to  regularise  the  problem  in  the  sense  of  Tikhonov. 

• It  is  often  convenient  to  sample  the  inverted  function  on  a grid  which  differs  from 
the  set  of  centres  Xi  of  the  basis.  With  a localised  Fourier  transform,  such  a sampling 
can  be  performed  without  significant  aliasing. 

The  asymptotic  estimate  4/to(.t)  = 0(l/||a;||m+(d+1V2)  may  be  derived  from  (3.1)  and 
estimates  Im  with  large  argument  (see  Eq.  A4,  [5]),  a fact  which  should  inform  our  choice 
of  m. 

4 Choice  of  parameters 

One  agreeable  feature  of  Lewitt’s  radial  functions  is  that  the  choice  of  parameters  of 
the  functions  correspond  in  a natural  way  to  the  balance  between  representation  quality 
and  efficiency  of  computation.  For  example,  the  asymptotic  rate  of  decay  of  the  Fourier 
transform  increases  with  m (see  above),  but  so  does  the  cost  of  the  calculation  of  Im. 
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A similar  choice  arises  when  the  centres  lie  on  a uniform  (square  or  triangular)  grid 
F.  Let  A denote  the  grid  spacing  of  such  a grid,  i.  e.,  the  minimum  distance  between 
distinct  centres  in  T.  It  is  desirable  that  the  grid  ratio  a/ A be  small,  as  this  results  in 
sparsity  of  the  discretisation.  As  a guide  to  fixing  the  values  of  a and  the  grid  ratio, 
Lewitt  suggests  the  error  in  quasi  interpolation  to  a constant,  the  error  with  which  the 
function 

9(x)  :=5^*(a0 

ier 

approximates  the  function  whose  constant  value  is  that  of  g at  the  centres  (edge  effects 
are  ignored  here).  In  Figure  2 the  root  mean  square  of  this  representation  error  (estimated 
numerically)  is  shown  for  the  square  planar  grid,  m = 2 and  a range  of  values  of  a and 
a/ A.  The  distinctive  “trenches”  in  the  error  can  be  explained  with  Poisson  summation 
formula  (see  [11]), 

^ ^(x  + An)  = “ ^ exp(27rf n.a;/A)^(7z/A).  (4.1) 

ne  Z2  nez2 

The  summand  for  n = 0 in  the  second  sum  is  'I'(O),  so  the  representation  error  depends 
only  on  the  values  of  T on  the  dual  grid,  Z2/A.  Provided  that  T decays  rapidly,  we 
would  expect  a small  error  when  T is  zero,  or  close  to  zero,  for  the  dual  grid-nodes  close 
to  the  origin. 

By  (3.1),  ’F(x)  is  zero  exactly  when 

Jm+d/2(V(2™\\x\\)2  ~ai)  = °> 

i.  e.,  for  radial  values  ||x||  = Rk, 

Rk  = (fc  = 1, 2, . . .), 

where  rjk  is  the  fc-th  zero  of  Jm+d/2-  Thus,  the  requirement  that  that  the  fc-th  zero  of 
#(x)  occurs  at  the  radius  of  the  closest  non- zero  dual  grid  node  (i.e.,  Rk  = 1/A)  is  a 
constraint  on  the  values  of  a and  a/ A 

a = \J  (27ra/  A)2  — rj2.  (4.2) 

The  contours  (4.2)  agree  well  with  the  trenches  evident  in  Figure  2.  With  the  same  intent 
we  can  require  that  the  Z-th  zero  of  T (x)  occur  at  the  radius  of  the  second  closest  dual 
grid  node  ( Ri  = \/2/A).  Points  satisfying  both  of  these  constraints  can  be  expected  to 
have  a particularly  small  representation  error.  In  Figure  2 these  favourable  choices  are 
labelled  k:l. 

The  above  argument  can  be  also  be  applied  the  triangular  grid.  Establishing  the  Pois- 
son summation  formula  for  such  is  straightforward  (either  generalised  from  VII  Section 
2 of  [11]  or  specialised  from  the  formula  for  topological  groups  in  [6]),  and  one  finds  that 
dual  grid  is  the  triangular  grid  with  node  spacing  2/(A\/3).  The  representation  error  is, 
qualitatively,  similar  to  that  shown  in  Figure  2.  To  make  a quantitative  comparison  we 
plot,  in  Figure  3,  the  representation  error  on  the  principal  trench  (i.  e.,  along  the  contour 
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Fig.  2.  The  representation  error  of  the  square  planar  grid  for  m = 2.  The  lower  contour 
map  shows  the  root  mean  square  error  in  representation  for  different  values  of  the  grid 
ratio  a/  A and  localisation  a.  The  upper  figure  shows  the  error  along  the  trenches  evident 
in  the  lower.  Favourable  choices  of  the  parameters  are  marked  1:2.  1:3,...,  and  are  also 
shown  in  the  lower  figure.  Values  of  the  parameters  to  the  left  of  the  dashed  line  give 
rise  to  a diagonally  dominant  interpolation  matrix. 
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Fig.  3.  Error  in  the  principal  trench  for  square  (dashed)  and  triangular  (solid)  grids. 

(4.2)  for  k = 1 in  the  case  of  the  square  grid)  for  each  grid  type  and  a number  of  values 
of  m. 

To  ensure  a fair  comparison,  the  horizontal  scale  in  Figure  3 is  adjusted  for  each 
grid-type  to  give  equal  node  densities.  As  is  seen,  the  two  grid-types  have  similar  error 
performance,  suggesting  that  the  square  grid  (with  attendant  ease  of  implementation) 
is  to  be  preferred  in  practice. 

5 The  functions  of  Wendland 

It  is  interesting  to  compare  Lewitt ’s  functions  with  the  radial  basis  functions  of  Wendland 
[1,  14],  positive  definite  functions  whose  window  functions  are  piecewise  polynomial.  The 
positive  definiteness  of  Wendland’s  functions  indicate  their  usefulness  in  approximation, 
for  which  extensive  results  exist,  and  a number  of  recent  papers  have  explored  their  use 
in  the  discretisation  of  partial  differential  equations. 

The  use  of  Wendland’s  functions  in  x-ray  problems  does  not  appear  to  have  been 
investigated,  although  their  Abel  transforms  can  be  obtained  analytically.  We  do  not  ad- 
dress this  question  here,  but  indicate  why  Lewitt’s  functions  may  offer  some  advantages 
for  such  problems.  The  Fourier  transform  of  Wendland’s  function  $2,0,  whose  window  is 
<p2,o (r)  = (1  — r)+,  is  proportional  to 

2tt  r 

(27rr  — t)2tJ0(t)  dt  = 0(r~3)  (r  = ||a:||) 

(see  Section  3,  [14]).  In  Figure  4,  4>2,o  is  plotted  along  with  the  Fourier  transform  T2,  of 
Lewitt’s  function  with  a = 1 and  the  parameter  choice  1 :2  of  Figure  2.  Although  both 
have  the  same  asymptotic  decay  of  the  Fourier  transform,  Lewitt’s  is  more  localised 
about  zero  and  thus  may  offer  better  suppression  of  ghosts  in  x-ray  problems. 

Finally  we  mention  that  Buhmann  has  shown,  in  [1],  that  Wendland’s  window  func- 


J.  J.  Green 


Fig.  4.  Fourier  transforms  of  basis  functions, 
tion  admits  a convolution  representation  of  the  form 


roc 

p{r):=  {l-r2/t)\tng{t)dt  (5.1) 

Jo 

for  the  weight  g(t ) = (1  - 1)+  with  suitable  k and  n.  We  note  that  (5.1)  may  be  solved 
for  g,  since  substituting  x = r2  in  (5.1)  allows  it  to  be  reduced  to  a standard  integral 
equation  whose  solution, 

9 (*)  = f(x)  = p(r), 


can  be  found  in  Article  1.1-4.32  of  [10].  In  the  case  that  p is  Lewitt’s  window  ipm,  one 
may  use  the  differentiation  formula,  All  of  [5],  to  find  the  corresponding  weight  g.  For 
n = 1 we  find  that 


P(*)  = 


1 lm— 1(^) 
2 am~2  Im(a) 


(r), 


a weight  qualitatively  different  from  that  of  Wendland’s  function. 
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